Top Modules Exploratory Analysis

#setwd('/afs/inf.ed.ac.uk/user/s17/s1725186/Documents/PhD-Models/FirstPUModel/RMarkdowns')

library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(dendextend)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally)
library(expss)
library(polycor)
library(foreach) ; library(doParallel)
library(knitr)
library(biomaRt)
library(anRichment) ; library(BrainDiseaseCollection)
suppressWarnings(suppressMessages(library(WGCNA)))

SFARI_colour_hue = function(r) {
  pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}

Load preprocessed dataset (preprocessing code in 19_10_14_data_preprocessing.Rmd) and clustering (pipeline in 19_10_21_WGCNA.Rmd)

# Gandal dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame


# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)


# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_08-29-2019_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]


# Clusterings
clusterings = read_csv('./../Data/clusters.csv')


# Update DE_info with SFARI and Neuronal information
genes_info = DE_info %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>% 
  mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
  left_join(GO_neuronal, by='ID') %>% left_join(clusterings, by='ID') %>%
  mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
  mutate(gene.score=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`), 
         significant=padj<0.05 & !is.na(padj))

# Add gene symbol
getinfo = c('ensembl_gene_id','external_gene_id')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl',
               host='feb2014.archive.ensembl.org') ## Gencode v19
gene_names = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=genes_info$ID, mart=mart)

genes_info = genes_info %>% left_join(gene_names, by=c('ID'='ensembl_gene_id'))


clustering_selected = 'DynamicHybrid'
genes_info$Module = genes_info[,clustering_selected]

dataset = read.csv(paste0('./../Data/dataset_', clustering_selected, '.csv'))
dataset$Module = dataset[,clustering_selected]


rm(DE_info, GO_annotations, clusterings, getinfo, mart, dds)

Relation to external clinical traits

Quantifying module-trait associations

Using the hetcor function, that calculates Pearson, polyserial or polychoric correlations depending on the type of variables involved.

datTraits = datMeta %>% dplyr::select(Diagnosis, Region, Sex, Age, PMI, RNAExtractionBatch) %>%
            dplyr::rename('ExtractionBatch' = RNAExtractionBatch)

# Recalculate MEs with color labels
ME_object = datExpr %>% t %>% moduleEigengenes(colors = genes_info$Module)
MEs = orderMEs(ME_object$eigengenes)

# Calculate correlation between eigengenes and the traits and their p-values
moduleTraitCor = MEs %>% apply(2, function(x) hetcor(x, datTraits)$correlations[1,-1]) %>% t
rownames(moduleTraitCor) = colnames(MEs)
colnames(moduleTraitCor) = colnames(datTraits)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nrow(datExpr))

# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)

# In case there are any NAs
if(sum(!complete.cases(moduleTraitCor))>0){
  print(paste0(sum(is.na(moduleTraitCor)),' correlation(s) could not be calculated')) 
}

## [1] "1 correlation(s) could not be calculated"

rm(ME_object)

Note: The correlation between Module #F464E5 and Diagonsis is the one that cannot be calculated, weirdly enough, the thing that causes the error is that the initial correlation is too high, so it would be a very bad thing to lose this module because of this numerical error. I’m going to fill in its value using the polyserial function, which doesn’t give exactly the same results as the hetcor() function, but it’s quite similar.

# Calculate the correlation tha failed with hetcor()
moduleTraitCor['ME#F464E5','Diagnosis'] = polyserial(MEs[,'ME#F464E5'], datTraits$Diagnosis)

## Warning in polyserial(MEs[, "ME#F464E5"], datTraits$Diagnosis): initial
## correlation inadmissible, 1.01199040044846, set to 0.9999

I’m going to select all the modules that have an absolute correlation higher than 0.9 with Diagnosis to study them

# Sort moduleTraitCor by Diagnosis
moduleTraitCor = moduleTraitCor[order(moduleTraitCor[,1], decreasing=TRUE),]
moduleTraitPvalue = moduleTraitPvalue[order(moduleTraitCor[,1], decreasing=TRUE),]

# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)


labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels =  gsub('ME','',rownames(moduleTraitCor)), 
               yColorWidth=0, colors = brewer.pal(11,'PiYG'), bg.lab.y = gsub('ME','',rownames(moduleTraitCor)),
               textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.8, cex.lab.y = 0.75, zlim = c(-1,1),
               main = paste('Module-Trait relationships'))

diagnosis_cor = data.frame('Module' = gsub('ME','',rownames(moduleTraitCor)),
                           'MTcor' = moduleTraitCor[,'Diagnosis'],
                           'MTpval' = moduleTraitPvalue[,'Diagnosis'])

genes_info = genes_info %>% left_join(diagnosis_cor, by='Module')

rm(moduleTraitPvalue, datTraits, textMatrix, diagnosis_cor)

Studying the modules with the highest absolute correlation to Diagnosis

The modules consist mainly of points with very high (absolute) values in PC2 (which we know is related to lfc), so this result is consistent with the high correlation between Module and Diagnosis, although some of the points with the highest PC2 values do not belong to these top modules

top_modules = gsub('ME','',rownames(moduleTraitCor)[abs(moduleTraitCor[,'Diagnosis'])>0.9])

cat(paste0('Top modules selected: ', paste(top_modules, collapse=', '),'\n'))

## Top modules selected: #F464E5, #00BECA

pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
            dplyr::select(ID, external_gene_id, PC1, PC2, Module, gene.score) %>%
            mutate(ImportantModules = ifelse(Module %in% top_modules, as.character(Module), 'Others')) %>%
            mutate(color = ifelse(ImportantModules=='Others','gray',ImportantModules),
                   alpha = ifelse(ImportantModules=='Others', 0.2, 0.4),
                   gene_id = paste0(ID, ' (', external_gene_id, ')'))

table(plot_data$ImportantModules)

## 
## #00BECA #F464E5  Others 
##    1638     626   13884

ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=ImportantModules)) + 
         geom_point(alpha=plot_data$alpha, color=plot_data$color, aes(ID=gene_id)) + theme_minimal() + 
           ggtitle('Modules with strongest relation to Diagnosis'))

rm(pca)

---

Module Membership vs Gene Significance

create_plot = function(module){
  
  plot_data = dataset %>% dplyr::select(ID, paste0('MM.',gsub('#','',module)), GS, gene.score) %>% filter(dataset$Module==module)
  colnames(plot_data)[2] = 'Module'
  
  SFARI_colors = as.numeric(names(table(as.character(plot_data$gene.score)[plot_data$gene.score!='None'])))
  
  p = ggplotly(plot_data %>% ggplot(aes(Module, GS, color=gene.score)) + geom_point(alpha=0.5, aes(ID=ID)) + ylab('Gene Significance') +
               scale_color_manual(values=SFARI_colour_hue(r=c(SFARI_colors,8))) + theme_minimal() + xlab('Module Membership') +
               ggtitle(paste0('Module ', module,'  (MTcor = ', round(moduleTraitCor[paste0('ME',module),1],2),')')))
  
  return(p)
}

create_plot(top_modules[1])

create_plot(top_modules[2])

rm(create_plot)

---

SFARI Genes

List of SFARI Genes in top modules ordered by SFARI score and Gene Significance

table_data = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
             dplyr::select(ID, external_gene_id, GS, gene.score, Module) %>% arrange(gene.score, desc(abs(GS))) %>%
             dplyr::rename('Ensembl ID'=ID, 'Gene Symbol'=external_gene_id, 
                    'SFARI score'=gene.score, 'Gene Significance'=GS)


kable(table_data %>% filter(Module == top_modules[1] & `SFARI score` != 'None') %>% dplyr::select(-Module),
      caption=paste0('SFARI Genes for Module ', top_modules[1]))

SFARI Genes for Module #F464E5

Ensembl ID	Gene Symbol	Gene Significance	SFARI score
ENSG00000038382	TRIO	0.6264110	2
ENSG00000168769	TET2	0.8178650	3
ENSG00000083168	KAT6A	0.5998865	3
ENSG00000113742	CPEB4	0.3701387	3
ENSG00000145020	AMT	0.3424890	3
ENSG00000132510	KDM6B	0.3330706	3
ENSG00000128573	FOXP2	0.2489053	3
ENSG00000112902	SEMA5A	0.2191666	3
ENSG00000165995	CACNB2	0.1701529	3
ENSG00000137801	THBS1	0.8076180	4
ENSG00000106366	SERPINE1	0.7907461	4
ENSG00000011422	PLAUR	0.7753875	4
ENSG00000172554	SNTG2	0.7354320	4
ENSG00000181481	RNF135	0.6768553	4
ENSG00000152208	GRID2	0.6205857	4
ENSG00000227184	EPPK1	0.5946476	4
ENSG00000183098	GPC6	0.5839280	4
ENSG00000188785	ZNF548	0.5801457	4
ENSG00000198589	LRBA	0.5579943	4
ENSG00000182372	CLN8	0.3972145	4
ENSG00000071246	VASH1	0.3828434	4
ENSG00000196277	GRM7	0.3773957	4
ENSG00000239389	PCDHA13	0.3419447	4
ENSG00000007314	SCN4A	0.3327403	4
ENSG00000182472	CAPN12	0.3178125	4
ENSG00000224389	C4B	0.3077509	4
ENSG00000215045	GRID2IP	0.2888427	4
ENSG00000130508	PXDN	0.1235447	4
ENSG00000187957	DNER	0.0565791	4
ENSG00000100024	UPB1	0.3133084	5
ENSG00000026508	CD44	0.2637995	5
ENSG00000136244	IL6	0.1960508	5
ENSG00000184588	PDE4B	0.1947772	5
ENSG00000141526	SLC16A3	0.8254828	6
ENSG00000072364	AFF4	0.7711840	6
ENSG00000171791	BCL2	0.3611981	6

kable(table_data %>% filter(Module == top_modules[2] & `SFARI score` != 'None') %>% dplyr::select(-Module),
      caption=paste0('SFARI Genes for Module ', top_modules[2]))

SFARI Genes for Module #00BECA

Ensembl ID	Gene Symbol	Gene Significance	SFARI score
ENSG00000136535	TBR1	-0.6544304	1
ENSG00000136531	SCN2A	-0.5379355	1
ENSG00000174469	CNTNAP2	-0.7002325	2
ENSG00000061676	NCKAP1	-0.6918489	2
ENSG00000139613	SMARCC2	-0.4332840	2
ENSG00000196557	CACNA1H	-0.4139346	2
ENSG00000119866	BCL11A	-0.3842245	2
ENSG00000157445	CACNA2D3	-0.3391182	2
ENSG00000144619	CNTN4	-0.3111072	2
ENSG00000074590	NUAK1	-0.8234159	3
ENSG00000144285	SCN1A	-0.8004590	3
ENSG00000196876	SCN8A	-0.7890809	3
ENSG00000170579	DLGAP1	-0.7792697	3
ENSG00000005955	GGNBP2	-0.7656571	3
ENSG00000132294	EFR3A	-0.7504244	3
ENSG00000078328	RBFOX1	-0.7315557	3
ENSG00000003147	ICA1	-0.7117475	3
ENSG00000197535	MYO5A	-0.7116705	3
ENSG00000175497	DPP10	-0.7061637	3
ENSG00000182621	PLCB1	-0.6921795	3
ENSG00000171759	PAH	-0.6897098	3
ENSG00000021645	NRXN3	-0.6713796	3
ENSG00000166501	PRKCB	-0.6155044	3
ENSG00000183454	GRIN2A	-0.5294268	3
ENSG00000146830	GIGYF1	-0.5256914	3
ENSG00000185345	PARK2	-0.4980108	3
ENSG00000164506	STXBP5	-0.4858781	3
ENSG00000132024	CC2D1A	-0.4128981	3
ENSG00000050030	KIAA2022	-0.4002974	3
ENSG00000168116	KIAA1586	-0.3735057	3
ENSG00000133026	MYH10	-0.3267293	3
ENSG00000101489	CELF4	-0.3181586	3
ENSG00000185008	ROBO2	-0.3077778	3
ENSG00000182256	GABRG3	-0.2873575	3
ENSG00000139174	PRICKLE1	-0.2147716	3
ENSG00000138411	HECW2	-0.1373931	3
ENSG00000140945	CDH13	0.0845059	3
ENSG00000166147	FBN1	-0.0824793	3
ENSG00000171723	GPHN	-0.0749364	3
ENSG00000149972	CNTN5	-0.0485789	3
ENSG00000144406	UNC80	-0.8620647	4
ENSG00000135631	RAB11FIP5	-0.8335775	4
ENSG00000132639	SNAP25	-0.8106671	4
ENSG00000163618	CADPS	-0.8086094	4
ENSG00000159082	SYNJ1	-0.7858980	4
ENSG00000104093	DMXL2	-0.7629031	4
ENSG00000172260	NEGR1	-0.7083293	4
ENSG00000163399	ATP1A1	-0.6896505	4
ENSG00000115840	SLC25A12	-0.6828765	4
ENSG00000196090	PTPRT	-0.6742079	4
ENSG00000130477	UNC13A	-0.6667346	4
ENSG00000154263	ABCA10	-0.6407178	4
ENSG00000153575	TUBGCP5	-0.6283954	4
ENSG00000136928	GABBR2	-0.6252262	4
ENSG00000144290	SLC4A10	-0.5769256	4
ENSG00000144331	ZNF385B	-0.5383310	4
ENSG00000186094	AGBL4	-0.5377404	4
ENSG00000107105	ELAVL2	-0.5374477	4
ENSG00000196433	ASMT	-0.5323119	4
ENSG00000155961	RAB39B	-0.5256637	4
ENSG00000116141	MARK1	-0.5062338	4
ENSG00000155886	SLC24A2	-0.4946020	4
ENSG00000109158	GABRA4	-0.4943945	4
ENSG00000198010	DLGAP2	-0.4942182	4
ENSG00000243232	PCDHAC2	-0.4777463	4
ENSG00000128594	LRRC4	-0.4719222	4
ENSG00000175161	CADM2	-0.4567470	4
ENSG00000182901	RGS7	-0.4466930	4
ENSG00000133019	CHRM3	-0.4273514	4
ENSG00000172915	NBEA	-0.4149285	4
ENSG00000156113	KCNMA1	-0.4061200	4
ENSG00000081189	MEF2C	-0.4000890	4
ENSG00000114374	USP9Y	-0.3520447	4
ENSG00000174442	ZWILCH	-0.3276604	4
ENSG00000136653	RASSF5	-0.3267914	4
ENSG00000151790	TDO2	-0.3195670	4
ENSG00000152495	CAMK4	-0.3178490	4
ENSG00000060140	STYK1	-0.3136339	4
ENSG00000064687	ABCA7	-0.3112037	4
ENSG00000165355	FBXO33	-0.3050639	4
ENSG00000135127	CCDC64	-0.3047331	4
ENSG00000170289	CNGB3	-0.2731444	4
ENSG00000139915	MDGA2	-0.2430880	4
ENSG00000152402	GUCY1A2	-0.2327724	4
ENSG00000111615	KRR1	-0.2248392	4
ENSG00000100359	SGSM3	-0.2075939	4
ENSG00000159640	ACE	-0.1924407	4
ENSG00000004468	CD38	-0.1674841	4
ENSG00000128513	POT1	-0.1434731	4
ENSG00000244588	RAD21L1	-0.1256607	4
ENSG00000165379	LRFN5	-0.1218335	4
ENSG00000100038	TOP3B	-0.0993690	4
ENSG00000152214	RIT2	0.0792829	4
ENSG00000165246	NLGN4Y	-0.0675768	4
ENSG00000181036	FCRL6	0.0604928	4
ENSG00000007168	PAFAH1B1	-0.8564947	5
ENSG00000087470	DNM1L	-0.8491182	5
ENSG00000104725	NEFL	-0.7912259	5
ENSG00000022355	GABRA1	-0.7777200	5
ENSG00000118596	SLC16A7	-0.7640017	5
ENSG00000171735	CAMTA1	-0.7302306	5
ENSG00000146469	VIP	-0.6852620	5
ENSG00000152969	JAKMIP1	-0.6814938	5
ENSG00000106123	EPHB6	-0.6655863	5
ENSG00000173726	TOMM20	-0.6303781	5
ENSG00000100362	PVALB	-0.4940183	5
ENSG00000007171	NOS2	-0.4656878	5
ENSG00000102468	HTR2A	-0.4562790	5
ENSG00000145979	TBC1D7	-0.4227883	5
ENSG00000130707	ASS1	-0.3913746	5
ENSG00000158258	CLSTN2	-0.3471046	5
ENSG00000196730	DAPK1	-0.2855914	5
ENSG00000166603	MC4R	-0.2783239	5
ENSG00000185666	SYN3	-0.2690820	5
ENSG00000157168	NRG1	-0.1662121	5
ENSG00000166866	MYO1A	-0.1473913	5
ENSG00000106113	CRHR2	-0.0441381	5
ENSG00000077279	DCX	-0.5750475	6
ENSG00000148680	HTR7	-0.4261417	6
ENSG00000125780	TGM3	-0.3834141	6

Modules with the strongest module-diagnosis correlation should have the highest percentage of SFARI Genes, but this doesn’t seem to be the case (even the opposite may be true)

plot_data = dataset %>% mutate('hasSFARIscore' = gene.score!='None') %>% 
            group_by(Module, MTcor, hasSFARIscore) %>% summarise(p=n()) %>% 
            left_join(dataset %>% group_by(Module) %>% summarise(n=n()), by='Module') %>% 
            mutate(p=round(p/n*100,2)) 

for(i in 1:nrow(plot_data)){
  this_row = plot_data[i,]
  if(this_row$hasSFARIscore==FALSE & this_row$p==100){
    new_row = this_row
    new_row$hasSFARIscore = TRUE
    new_row$p = 0
    plot_data = plot_data %>% rbind(new_row)
  }
}

plot_data = plot_data %>% filter(hasSFARIscore==TRUE)

ggplotly(plot_data %>% ggplot(aes(MTcor, p, size=n)) + geom_smooth(color='gray', se=FALSE) +
         geom_point(color=plot_data$Module, alpha=0.5, aes(id=Module)) + geom_hline(yintercept=mean(plot_data$p), color='gray') +
         xlab('Module-Diagnosis correlation') + ylab('% of SFARI genes') +
         theme_minimal() + theme(legend.position = 'none'))

rm(i, this_row, new_row, plot_data)

---

Module Eigengenes

Since these modules have the strongest relation to autism, this pattern should be reflected in their model eigengenes, having two different behaviours for the samples corresponding to autism and the ones corresponding to control.

In both cases, the Eigengenes separate the behaviour between autism and control samples very clearly!

plot_EGs = function(module){

  plot_data = data.frame('ID' = rownames(MEs), 'MEs' = MEs[,paste0('ME',module)], 'Diagnosis' = datMeta$Diagnosis)

  p = plot_data %>% ggplot(aes(Diagnosis, MEs, fill=Diagnosis)) + geom_boxplot() + theme_minimal() + theme(legend.position='none') +
                    ggtitle(paste0('Module ', module, '  (MTcor=',round(moduleTraitCor[paste0('ME',module),1],2),')'))
  return(p)
}

p1 = plot_EGs(top_modules[1])
p2 = plot_EGs(top_modules[2])

grid.arrange(p1, p2, nrow=1)

rm(plot_EGs, p1, p2)

---

Identifying important genes

Selecting the modules with the highest correlation to Diagnosis, and, from them, the genes with the highest module membership-(absolute) gene significance

*Ordered by \(\frac{MM+|GS|}{2}\)

There aren’t that many SFARI genes in the top genes of the modules and not a single one belonging to scores 1 and 2

create_table = function(module){
  top_genes = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>% 
              dplyr::select(ID, external_gene_id, paste0('MM.',gsub('#','',module)), GS, gene.score) %>%
              filter(dataset$Module==module) %>% dplyr::rename('MM' = paste0('MM.',gsub('#','',module))) %>% 
              mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>% top_n(20)
  return(top_genes)
}

top_genes_1 = create_table(top_modules[1])
kable(top_genes_1, caption=paste0('Top 10 genes for module ', top_modules[1], '  (MTcor = ',
                                  round(moduleTraitCor[paste0('ME',top_modules[1]),1],2),')'))

Top 10 genes for module #F464E5 (MTcor = 1)

ID	external_gene_id	MM	GS	gene.score	importance
ENSG00000143384	MCL1	0.8476929	0.9075109	None	0.8776019
ENSG00000161638	ITGA5	0.8510043	0.8600246	None	0.8555144
ENSG00000089159	PXN	0.8085675	0.8694821	None	0.8390248
ENSG00000196935	SRGAP1	0.7501468	0.9223277	None	0.8362373
ENSG00000003402	CFLAR	0.7941123	0.8169199	None	0.8055161
ENSG00000148841	ITPRIP	0.7881631	0.8162584	None	0.8022107
ENSG00000124782	RREB1	0.7589325	0.8451372	None	0.8020348
ENSG00000106366	SERPINE1	0.8065142	0.7907461	4	0.7986301
ENSG00000150457	LATS2	0.7443403	0.8501398	None	0.7972400
ENSG00000138119	MYOF	0.7514562	0.8343142	None	0.7928852
ENSG00000073792	IGF2BP2	0.7392534	0.8389676	None	0.7891105
ENSG00000138166	DUSP5	0.8112245	0.7636771	None	0.7874508
ENSG00000072364	AFF4	0.7963444	0.7711840	6	0.7837642
ENSG00000120278	PLEKHG1	0.7556018	0.8087123	None	0.7821570
ENSG00000162745	OLFML2B	0.7023713	0.8526646	None	0.7775179
ENSG00000154640	BTG3	0.7244112	0.8301676	None	0.7772894
ENSG00000133639	BTG1	0.7423866	0.8079258	None	0.7751562
ENSG00000168769	TET2	0.7312714	0.8178650	3	0.7745682
ENSG00000120690	ELF1	0.7613055	0.7876526	None	0.7744790
ENSG00000158615	PPP1R15B	0.7289451	0.8197298	None	0.7743375

top_genes_2 = create_table(top_modules[2])
kable(top_genes_2, caption=paste0('Top 10 genes for module ', top_modules[2], '  (MTcor = ',
                                  round(moduleTraitCor[paste0('ME',top_modules[2]),1],2),')'))

Top 10 genes for module #00BECA (MTcor = -0.92)

ID	external_gene_id	MM	GS	gene.score	importance
ENSG00000050748	MAPK9	0.9055266	-0.9401065	None	0.9228165
ENSG00000108395	TRIM37	0.9056872	-0.9297206	None	0.9177039
ENSG00000138078	PREPL	0.8926949	-0.9043456	None	0.8985203
ENSG00000177432	NAP1L5	0.8738603	-0.9048896	None	0.8893750
ENSG00000163577	EIF5A2	0.8622485	-0.9110526	None	0.8866505
ENSG00000128881	TTBK2	0.8938083	-0.8754731	None	0.8846407
ENSG00000176490	DIRAS1	0.8610032	-0.8968398	None	0.8789215
ENSG00000171132	PRKCE	0.8477771	-0.9090963	None	0.8784367
ENSG00000155097	ATP6V1C1	0.8504461	-0.8966833	None	0.8735647
ENSG00000196876	SCN8A	0.9470763	-0.7890809	3	0.8680786
ENSG00000114573	ATP6V1A	0.8139207	-0.9175213	None	0.8657210
ENSG00000111674	ENO2	0.8761932	-0.8485455	None	0.8623694
ENSG00000125814	NAPB	0.8922318	-0.8313011	None	0.8617664
ENSG00000184368	MAP7D2	0.8332534	-0.8807598	None	0.8570066
ENSG00000131437	KIF3A	0.8671455	-0.8464587	None	0.8568021
ENSG00000144285	SCN1A	0.9127631	-0.8004590	3	0.8566110
ENSG00000162694	EXTL2	0.8023344	-0.8984415	None	0.8503880
ENSG00000172348	RCAN2	0.7828069	-0.9152892	None	0.8490480
ENSG00000130540	SULT4A1	0.8830724	-0.8116827	None	0.8473775
ENSG00000132639	SNAP25	0.8784733	-0.8106671	4	0.8445702

rm(create_table)

pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(color = ifelse(Module %in% top_modules, as.character(Module), 'gray')) %>%
            mutate(alpha = ifelse(color %in% top_modules & 
                                  ID %in% c(as.character(top_genes_1$ID), 
                                            as.character(top_genes_2$ID)), 1, 0.1))

plot_data %>% ggplot(aes(PC1, PC2)) + geom_point(alpha=plot_data$alpha, color=plot_data$color) + 
              theme_minimal() + ggtitle('Important genes identified through WGCNA')

---

Enrichment Analysis

Using the package anRichment

• It was designed by Peter Langfelder explicitly to perform enrichmen analysis on WGCNA’s modules in brain-related experiments (mainly Huntington’s Disease)

• It has packages with brain annotations:

  • BrainDiseaseCollection: A Brain Disease Gene Set Collection for anRichment

  • MillerAIBSCollection: (included in anRichment) Contains gene sets collected by Jeremy A. Miller at AIBS of various cell type and brain region marker sets, gene sets collected from expression studies of developing brain, as well as a collection of transcription factor (TF) targets from the original ChEA study

• The tutorial says it’s an experimental package

• It’s not on CRAN nor in Bioconductor

# Prepare dataset

# Create dataset with top modules membership and removing the genes without an assigned module
EA_dataset = data.frame('ensembl_gene_id' = genes_info$ID,
                        module = ifelse(genes_info$Module %in% top_modules, genes_info$Module, 'other')) %>%
             filter(genes_info$Module!='gray')

# Assign Entrez Gene Id to each gene
getinfo = c('ensembl_gene_id','entrezgene')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', host='feb2014.archive.ensembl.org')
biomart_output = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=EA_dataset$ensembl_gene_id, mart=mart)

## Cache found

EA_dataset = EA_dataset %>% left_join(biomart_output, by='ensembl_gene_id')

for(tm in top_modules){
  cat(paste0('\n',sum(EA_dataset$module==tm & is.na(EA_dataset$entrezgene)), ' genes from top module ',
             tm, ' don\'t have an Entrez Gene ID')) 
}

## 
## 8 genes from top module #F464E5 don't have an Entrez Gene ID
## 29 genes from top module #00BECA don't have an Entrez Gene ID

rm(getinfo, mart, biomart_output, tm)

# Manual: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/GeneAnnotation/Tutorials/anRichment-Tutorial1.pdf

collectGarbage()

# EA_dataset = rbind(EA_dataset[EA_dataset$module!='other',], EA_dataset[EA_dataset$module=='other',][sample(sum(EA_dataset$module=='other'), 1000),])

# Prepare datasets
GO_col = buildGOcollection(organism = 'human', verbose = 0)

## Loading required package: org.Hs.eg.db

## 

## 'select()' returned 1:many mapping between keys and columns

## 'select()' returned 1:1 mapping between keys and columns

internal_col = internalCollection(organism = 'human')
MillerAIBS_col = MillerAIBSCollection(organism = 'human')
BrainDisease_col = BrainDiseaseCollection(organism = 'human')
combined_col = mergeCollections(GO_col, internal_col, MillerAIBS_col, BrainDisease_col)

# Print collections used
cat('Using collections: ')

## Using collections:

knownGroups(combined_col, sortBy = 'size')

##  [1] "GO"                                         
##  [2] "GO.BP"                                      
##  [3] "GO.MF"                                      
##  [4] "GO.CC"                                      
##  [5] "JA Miller at AIBS"                          
##  [6] "Chip-X enrichment analysis (ChEA)"          
##  [7] "Brain"                                      
##  [8] "JAM"                                        
##  [9] "Prenatal brain"                             
## [10] "Brain region markers"                       
## [11] "Cortex"                                     
## [12] "Brain region marker enriched gene sets"     
## [13] "WGCNA"                                      
## [14] "BrainRegionMarkers"                         
## [15] "BrainRegionMarkers.HBA"                     
## [16] "BrainRegionMarkers.HBA.localMarker(top200)" 
## [17] "Postnatal brain"                            
## [18] "ImmunePathways"                             
## [19] "Markers of cortex layers"                   
## [20] "BrainLists"                                 
## [21] "Cell type markers"                          
## [22] "Germinal brain"                             
## [23] "BrainRegionMarkers.HBA.globalMarker(top200)"
## [24] "Accelerated evolution"                      
## [25] "Postmitotic brain"                          
## [26] "BrainLists.Blalock_AD"                      
## [27] "BrainLists.DiseaseGenes"                    
## [28] "BloodAtlases"                               
## [29] "Verge Disease Genes"                        
## [30] "BloodAtlases.Whitney"                       
## [31] "BrainLists.JAXdiseaseGene"                  
## [32] "BrainLists.MO"                              
## [33] "Age-associated genes"                       
## [34] "BrainLists.Lu_Aging"                        
## [35] "Cell type marker enriched gene sets"        
## [36] "BrainLists.CA1vsCA3"                        
## [37] "BrainLists.MitochondrialType"               
## [38] "BrainLists.MO.2+_26Mar08"                   
## [39] "BrainLists.MO.Sugino"                       
## [40] "BloodAtlases.Gnatenko2"                     
## [41] "BloodAtlases.Kabanova"                      
## [42] "BrainLists.Voineagu"                        
## [43] "StemCellLists"                              
## [44] "StemCellLists.Lee"

# Perform Enrichment Analysis
enrichment = enrichmentAnalysis(classLabels = EA_dataset$module, identifiers = EA_dataset$entrezgene,
                                refCollection = combined_col, useBackground = 'given', 
                                threshold = 1e-4, thresholdType = 'Bonferroni',
                                getOverlapEntrez = FALSE, getOverlapSymbols = TRUE)

##  enrichmentAnalysis: preparing data..
##   ..working on label set 1 ..

Results

kable(enrichment$enrichmentTable %>% filter(class==top_modules[1]) %>% 
      dplyr::select(dataSetID, shortDataSetName, inGroups, Bonferroni, FDR,
                    effectiveClassSize, effectiveSetSize, nCommonGenes),
      caption = paste0('Enriched terms for module ', top_modules[1], ' (MTcor = ',
                       round(genes_info$MTcor[genes_info$Module==top_modules[1]][1],4), ')'))

Enriched terms for module #F464E5 (MTcor = 0.9999)

dataSetID	shortDataSetName	inGroups	Bonferroni	FDR	effectiveClassSize	effectiveSetSize	nCommonGenes
JAMiller.AIBS.000434	Genes bound by RELA in HUMAN FIBROSARCOMA from PMID 24523406	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	0.0000002	0.0e+00	642	1043	89
GO:0001525	angiogenesis	GO|GO.BP	0.0000071	1.0e-07	642	434	48
GO:0072358	cardiovascular system development	GO|GO.BP	0.0000273	3.0e-07	642	624	59
GO:0007166	cell surface receptor signaling pathway	GO|GO.BP	0.0000365	3.0e-07	642	2498	157
GO:0001944	vasculature development	GO|GO.BP	0.0000429	4.0e-07	642	615	58
GO:0072359	circulatory system development	GO|GO.BP	0.0000488	4.0e-07	642	941	77
GO:0001568	blood vessel development	GO|GO.BP	0.0001549	1.3e-06	642	587	55
GO:0048514	blood vessel morphogenesis	GO|GO.BP	0.0002049	1.6e-06	642	511	50
GO:0035295	tube development	GO|GO.BP	0.0005272	3.7e-06	642	884	71
GO:0030198	extracellular matrix organization	GO|GO.BP	0.0007586	5.1e-06	642	316	36

kable(enrichment$enrichmentTable %>% filter(class==top_modules[2]) %>% 
      dplyr::select(dataSetID, shortDataSetName, inGroups, Bonferroni, FDR,
                    effectiveClassSize, effectiveSetSize, nCommonGenes),
      caption = paste0('Enriched terms for module ', top_modules[2], ' (MTcor = ',
                       round(genes_info$MTcor[genes_info$Module==top_modules[2]][1],4), ')'))

Enriched terms for module #00BECA (MTcor = -0.9186)

dataSetID	shortDataSetName	inGroups	Bonferroni	FDR	effectiveClassSize	effectiveSetSize	nCommonGenes
JAMiller.AIBS.000142	Highest in CP of 13-16 post-conception weeks human	JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex|Postmitotic brain	0.00e+00	0e+00	1637	1218	301
JAMiller.AIBS.000052	CortexWGCNA 15-21 post-conception weeks C26	JA Miller at AIBS|Brain|Prenatal brain|Cortex|WGCNA	0.00e+00	0e+00	1637	726	189
JAMiller.AIBS.000150	Highest in CP of E14.5 mouse	JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex|Postmitotic brain	0.00e+00	0e+00	1637	1276	273
JAMiller.AIBS.000155	Lowest in VZ of E14.5 mouse	JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex	0.00e+00	0e+00	1637	1674	317
JAMiller.AIBS.000569	WGCNA humanSpecificOlivedrab2Module frontalCtx FOXP2	JA Miller at AIBS|Brain|Postnatal brain|Cortex|WGCNA	0.00e+00	0e+00	1637	4057	599
JAM:002744	Autism_differential_expression_across_at_least_one_comparison	JAM|BrainLists|BrainLists.Voineagu	0.00e+00	0e+00	1637	765	175
JAM:003016	downAD_synapticTransmission	JAM|BrainLists|BrainLists.Blalock_AD	0.00e+00	0e+00	1637	89	48
JAMiller.AIBS.000141	CP enriched in E14.5 mouse	JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex|Postmitotic brain	0.00e+00	0e+00	1637	569	139
JAMiller.AIBS.000506	Genes bound by SUZ12 in MOUSE MESC from PMID 20075857	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	0.00e+00	0e+00	1637	3380	497
GO:0097458	neuron part	GO|GO.CC	0.00e+00	0e+00	1637	1607	281
JAMiller.AIBS.000570	WGCNA Olivedrab2ModuleGenes with enriched ELAVL2 targets	JA Miller at AIBS|Brain|Postnatal brain|Cortex|WGCNA	0.00e+00	0e+00	1637	483	119
GO:0045202	synapse	GO|GO.CC	0.00e+00	0e+00	1637	1114	208
GO:0099536	synaptic signaling	GO|GO.BP	0.00e+00	0e+00	1637	652	142
GO:0044456	synapse part	GO|GO.CC	0.00e+00	0e+00	1637	893	175
GO:0007268	chemical synaptic transmission	GO|GO.BP	0.00e+00	0e+00	1637	638	138
GO:0098916	anterograde trans-synaptic signaling	GO|GO.BP	0.00e+00	0e+00	1637	638	138
GO:0099537	trans-synaptic signaling	GO|GO.BP	0.00e+00	0e+00	1637	646	139
JAMiller.AIBS.000005	CPi markers at 21 post-conception weeks	JA Miller at AIBS|Brain|Prenatal brain|Cortex|Markers of cortex layers|Postmitotic brain	0.00e+00	0e+00	1637	312	85
JAM:003072	Tail of Caudate Nucleus_IN_Striatum	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	162	57
JAM:002805	Cerebral Cortex	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.globalMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	161	56
JAMiller.AIBS.000149	Lowest in VZ of 13-16 post-conception weeks human	JA Miller at AIBS|Brain|Prenatal brain|Markers of cortex layers|Cortex	0.00e+00	0e+00	1637	618	131
JAMiller.AIBS.000123	HippocampusWGCNA turquoise DGenriched upAge	JA Miller at AIBS|Brain|Postnatal brain|WGCNA	0.00e+00	0e+00	1637	1104	198
JAM:002751	Basal Pons	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.globalMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	167	56
GO:0043005	neuron projection	GO|GO.CC	0.00e+00	0e+00	1637	1214	210
JAM:003054	subiculum_IN_Hippocampal Formation	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	164	53
JAM:002739	arcuate nucleus of medulla_IN_Myelencephalon	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	168	53
GO:0097060	synaptic membrane	GO|GO.CC	0.00e+00	0e+00	1637	412	94
JAM:002907	inferior olivary complex_IN_Myelencephalon	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	170	53
JAMiller.AIBS.000504	Genes bound by SUZ12 in mouse MESC from PMID 18692474	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	0.00e+00	0e+00	1637	1362	222
JAM:002824	Dentate Nucleus_IN_Cerebellar Nucleus	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	165	51
GO:0030424	axon	GO|GO.CC	0.00e+00	0e+00	1637	574	116
GO:1902495	transmembrane transporter complex	GO|GO.CC	0.00e+00	0e+00	1637	289	72
JAM:003058	Substantia Nigra, pars compacta_IN_Mesencephalon	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	169	51
GO:0034702	ion channel complex	GO|GO.CC	0.00e+00	0e+00	1637	269	68
JAMiller.AIBS.000364	Genes bound by MTF2 in MOUSE MESC from PMID 20144788	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	0.00e+00	0e+00	1637	2292	330
JAM:002918	lateral medullary reticular group_IN_Myelencephalon	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	0.00e+00	0e+00	1637	162	49
GO:1990351	transporter complex	GO|GO.CC	1.00e-07	0e+00	1637	297	72
JAM:002764	downAging_mitochondria_synapse	JAM|BrainLists|BrainLists.Lu_Aging	1.00e-07	0e+00	1637	391	86
GO:0045211	postsynaptic membrane	GO|GO.CC	1.00e-07	0e+00	1637	309	73
JAMiller.AIBS.000503	Genes bound by SUZ12 in mouse MESC from PMID 18555785	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	2.00e-07	0e+00	1637	765	138
JAM:002991	downAD_synapticTransmission	JAM|BrainLists|BrainLists.Blalock_AD	3.00e-07	0e+00	1637	97	35
JAM:002967	Occipital Lobe_IN_Cerebral Cortex	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	4.00e-07	0e+00	1637	149	45
GO:0050877	nervous system process	GO|GO.BP	5.00e-07	0e+00	1637	886	153
JAM:002964	Nucleus Accumbens_IN_Striatum	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	6.00e-07	0e+00	1637	161	47
GO:0098794	postsynapse	GO|GO.CC	6.00e-07	0e+00	1637	591	113
GO:0120025	plasma membrane bounded cell projection	GO|GO.CC	8.00e-07	0e+00	1637	1919	280
GO:0034703	cation channel complex	GO|GO.CC	9.00e-07	0e+00	1637	203	54
JAM:002882	Hippocampal Formation	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.globalMarker(top200)|Brain region markers|Brain region marker enriched gene sets	1.00e-06	0e+00	1637	169	48
JAM:002763	downAD_metalIonTransport_glycoprotein	JAM|BrainLists|BrainLists.Blalock_AD	1.10e-06	0e+00	1637	283	67
JAM:003085	trochlear nucleus_IN_Mesencephalon	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	1.20e-06	0e+00	1637	153	45
GO:0007267	cell-cell signaling	GO|GO.BP	1.60e-06	0e+00	1637	1439	221
GO:0042995	cell projection	GO|GO.CC	2.60e-06	0e+00	1637	1989	286
GO:0050804	modulation of chemical synaptic transmission	GO|GO.BP	4.40e-06	1e-07	1637	411	85
GO:0005886	plasma membrane	GO|GO.CC	4.50e-06	1e-07	1637	4277	541
GO:0099177	regulation of trans-synaptic signaling	GO|GO.BP	5.00e-06	1e-07	1637	412	85
GO:0098793	presynapse	GO|GO.CC	6.10e-06	1e-07	1637	462	92
GO:0071944	cell periphery	GO|GO.CC	7.10e-06	1e-07	1637	4381	551
GO:0098978	glutamatergic synapse	GO|GO.CC	7.70e-06	1e-07	1637	341	74
JAM:002920	Lateral Nucleus_IN_Amygdala	JAM|BrainRegionMarkers|BrainRegionMarkers.HBA|BrainRegionMarkers.HBA.localMarker(top200)|Brain region markers|Brain region marker enriched gene sets	8.40e-06	1e-07	1637	167	46
GO:0005216	ion channel activity	GO|GO.MF	8.80e-06	1e-07	1637	362	77
JAM:002769	downAD_mitochondrion	JAM|BrainLists|BrainLists.Blalock_AD	1.04e-05	1e-07	1637	265	62
JAMiller.AIBS.000463	Genes bound by SMAD4 in HUMAN A2780 from PMID 21799915	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	1.18e-05	1e-07	1637	2085	294
JAMiller.AIBS.000042	CortexWGCNA 15-21 post-conception weeks C16 SPenriched	JA Miller at AIBS|Brain|Prenatal brain|Cortex|WGCNA	1.21e-05	1e-07	1637	198	51
GO:0036477	somatodendritic compartment	GO|GO.CC	1.63e-05	2e-07	1637	795	136
GO:0044463	cell projection part	GO|GO.CC	1.91e-05	2e-07	1637	1345	205
GO:0120038	plasma membrane bounded cell projection part	GO|GO.CC	1.91e-05	2e-07	1637	1345	205
GO:0006812	cation transport	GO|GO.BP	1.93e-05	2e-07	1637	1007	163
GO:0005261	cation channel activity	GO|GO.MF	2.04e-05	2e-07	1637	282	64
GO:0022803	passive transmembrane transporter activity	GO|GO.MF	2.14e-05	2e-07	1637	389	80
GO:0022838	substrate-specific channel activity	GO|GO.MF	2.96e-05	3e-07	1637	371	77
GO:0003008	system process	GO|GO.BP	3.92e-05	4e-07	1637	1505	223
JAM:002866	GlutamatergicNeuronsInMouseCortex_Sugino	JAM|BrainLists|BrainLists.MO|BrainLists.MO.Sugino|Cell type marker enriched gene sets|Brain region marker enriched gene sets	4.29e-05	4e-07	1637	367	76
GO:0015267	channel activity	GO|GO.MF	4.53e-05	4e-07	1637	388	79
JAMiller.AIBS.000436	Genes bound by REST in MOUSE MESC from PMID 21632747	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	4.54e-05	4e-07	1637	1675	243
JAMiller.AIBS.000505	Genes bound by SUZ12 in MOUSE MESC from PMID 18974828	JA Miller at AIBS|Chip-X enrichment analysis (ChEA)	5.82e-05	5e-07	1637	1387	208
GO:0043269	regulation of ion transport	GO|GO.BP	5.94e-05	5e-07	1637	613	110

Save Enrichment Analysis results

save(enrichment, file='./../Data/enrichmentAnalysis.RData')

Session info

sessionInfo()

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Scientific Linux 7.6 (Nitrogen)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] org.Hs.eg.db_3.10.0                      
##  [2] BrainDiseaseCollection_1.00              
##  [3] anRichment_1.01-2                        
##  [4] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
##  [5] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2  
##  [6] GenomicFeatures_1.38.2                   
##  [7] GenomicRanges_1.38.0                     
##  [8] GenomeInfoDb_1.22.0                      
##  [9] anRichmentMethods_0.90-1                 
## [10] WGCNA_1.68                               
## [11] fastcluster_1.1.25                       
## [12] dynamicTreeCut_1.63-1                    
## [13] GO.db_3.10.0                             
## [14] AnnotationDbi_1.48.0                     
## [15] IRanges_2.20.2                           
## [16] S4Vectors_0.24.3                         
## [17] Biobase_2.46.0                           
## [18] BiocGenerics_0.32.0                      
## [19] biomaRt_2.42.0                           
## [20] knitr_1.24                               
## [21] doParallel_1.0.15                        
## [22] iterators_1.0.12                         
## [23] foreach_1.4.7                            
## [24] polycor_0.7-10                           
## [25] expss_0.10.1                             
## [26] GGally_1.4.0                             
## [27] gridExtra_2.3                            
## [28] viridis_0.5.1                            
## [29] viridisLite_0.3.0                        
## [30] RColorBrewer_1.1-2                       
## [31] dendextend_1.13.3                        
## [32] plotly_4.9.2                             
## [33] glue_1.3.1                               
## [34] reshape2_1.4.3                           
## [35] forcats_0.4.0                            
## [36] stringr_1.4.0                            
## [37] dplyr_0.8.3                              
## [38] purrr_0.3.3                              
## [39] readr_1.3.1                              
## [40] tidyr_1.0.2                              
## [41] tibble_2.1.3                             
## [42] ggplot2_3.2.1                            
## [43] tidyverse_1.3.0                          
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1                backports_1.1.5            
##   [3] Hmisc_4.2-0                 BiocFileCache_1.10.2       
##   [5] plyr_1.8.5                  lazyeval_0.2.2             
##   [7] splines_3.6.0               crosstalk_1.0.0            
##   [9] BiocParallel_1.20.1         robust_0.4-18.2            
##  [11] digest_0.6.24               htmltools_0.4.0            
##  [13] fansi_0.4.1                 magrittr_1.5               
##  [15] checkmate_1.9.4             memoise_1.1.0              
##  [17] fit.models_0.5-14           cluster_2.0.8              
##  [19] annotate_1.64.0             Biostrings_2.54.0          
##  [21] modelr_0.1.5                matrixStats_0.55.0         
##  [23] askpass_1.1                 prettyunits_1.0.2          
##  [25] colorspace_1.4-1            blob_1.2.1                 
##  [27] rvest_0.3.5                 rappdirs_0.3.1             
##  [29] rrcov_1.4-7                 haven_2.2.0                
##  [31] xfun_0.8                    crayon_1.3.4               
##  [33] RCurl_1.95-4.12             jsonlite_1.6               
##  [35] genefilter_1.68.0           impute_1.60.0              
##  [37] survival_2.44-1.1           gtable_0.3.0               
##  [39] zlibbioc_1.32.0             XVector_0.26.0             
##  [41] DelayedArray_0.12.2         DEoptimR_1.0-8             
##  [43] scales_1.1.0                mvtnorm_1.0-11             
##  [45] DBI_1.1.0                   Rcpp_1.0.3                 
##  [47] xtable_1.8-4                progress_1.2.2             
##  [49] htmlTable_1.13.1            foreign_0.8-71             
##  [51] bit_1.1-15.2                preprocessCore_1.48.0      
##  [53] Formula_1.2-3               htmlwidgets_1.5.1          
##  [55] httr_1.4.1                  ellipsis_0.3.0             
##  [57] acepack_1.4.1               farver_2.0.3               
##  [59] pkgconfig_2.0.3             reshape_0.8.8              
##  [61] XML_3.99-0.3                nnet_7.3-12                
##  [63] dbplyr_1.4.2                locfit_1.5-9.1             
##  [65] later_1.0.0                 labeling_0.3               
##  [67] tidyselect_0.2.5            rlang_0.4.4                
##  [69] munsell_0.5.0               cellranger_1.1.0           
##  [71] tools_3.6.0                 cli_2.0.1                  
##  [73] generics_0.0.2              RSQLite_2.2.0              
##  [75] broom_0.5.4                 fastmap_1.0.1              
##  [77] evaluate_0.14               yaml_2.2.0                 
##  [79] bit64_0.9-7                 fs_1.3.1                   
##  [81] robustbase_0.93-5           nlme_3.1-139               
##  [83] mime_0.9                    xml2_1.2.2                 
##  [85] compiler_3.6.0              rstudioapi_0.10            
##  [87] curl_4.3                    reprex_0.3.0               
##  [89] geneplotter_1.64.0          pcaPP_1.9-73               
##  [91] stringi_1.4.6               highr_0.8                  
##  [93] lattice_0.20-38             Matrix_1.2-17              
##  [95] vctrs_0.2.2                 pillar_1.4.3               
##  [97] lifecycle_0.1.0             data.table_1.12.8          
##  [99] bitops_1.0-6                httpuv_1.5.2               
## [101] rtracklayer_1.46.0          R6_2.4.1                   
## [103] latticeExtra_0.6-28         promises_1.1.0             
## [105] codetools_0.2-16            MASS_7.3-51.4              
## [107] assertthat_0.2.1            SummarizedExperiment_1.16.1
## [109] DESeq2_1.26.0               openssl_1.4.1              
## [111] withr_2.1.2                 GenomicAlignments_1.22.1   
## [113] Rsamtools_2.2.2             GenomeInfoDbData_1.2.2     
## [115] hms_0.5.3                   grid_3.6.0                 
## [117] rpart_4.1-15                rmarkdown_1.14             
## [119] Cairo_1.5-10                shiny_1.4.0                
## [121] lubridate_1.7.4             base64enc_0.1-3